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Abstract 

The rapid development of signal processing on graphs provides a new perspective 
for processing large-scale data associated with irregular domains. In many practical 
applications, it is necessary to handle massive data sets through complex networks, 
in which most nodes have limited computing power. Designing efficient distributed 
algorithms is critical for this task. This paper focuses on the distributed reconstruction 
of a time-varying bandlimited graph signal based on observations sampled at a subset of 
selected nodes. A distributed least square reconstruction (DLSR) algorithm is proposed 
to recover the unknown signal iteratively, by allowing neighboring nodes to communicate 
with one another and make fast updates. DLSR uses a decay scheme to annihilate 
the out-of-band energy occurring in the reconstruction process, which is inevitably 
caused by the transmission delay in distributed systems. Proof of convergence and 
error bounds for DLSR are provided in this paper, suggesting that the algorithm is able 
to track time-varying graph signals and perfectly reconstruct time-invariant signals. 
The DLSR algorithm is numerically experimented with synthetic data and real-world 
sensor network data, which verifies its ability in tracking slowly time-varying graph 
signals. 

Keywords: Signal processing on graph, graph signal, distributed algorithm, sam¬ 
pling and reconstruction, time-varying signal. 


1 Introduction 

1.1 Signal Processing on Graph 

During the past few years, the emerging field of signal processing on graphs (see 0 ® has 
attracted vast research interests from multiple disciplines. Consider an IV-vertex undirected 
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graph, denoted as G{V,£), where V is the set of vertices with |V| = N and £ is the set 
of edges. A graph signal f E is a vector that assigns each vertex a real number. 
Equivalently, the vector f is often regarded as a function / : V i —> 3?. 

Graph-based signal processing has been developed to analyze data or signal associated 
with irregular domains, e.g., large-scale networks. It finds wide applications in sensor net¬ 
works [31, image processing [2], recommendation systems [Fj, etc. Existing research topics 
on graph signal processing include: graph signal sampling |6,[7j, uncertainty principle |8|, 
graph filtering [ 9 ], spectral graph wavelet PPl , graph signal compression 1121, graph signal 


multiresolution |13|14| , parametric dictionary learning 15 , graph signal coarsening 16.17 
etc. 


1.2 Problem Description and Related Works 

In this work, we study the distributed reconstruction of smooth graph signals based on 
sample measurements obtained at representative nodes. Suppose that the unknown graph 
signal lies in the low-frequency subspace, our reconstruction problem is to recover its missing 
entries from its known data/signal values sampled from a representative set of nodes. 

Some theoretical results have been established for the sampling problem of bandlimited 
graph-based signals; see e.g., [I8|-[20j. The relation between the sample size necessary to 
obtain unique reconstruction and the cutoff frequency of bandlimited signal space has been 
studied. Similar to classical results on time-domain irregular sampling, the idea of “frame” 
has been introduced for graph signal processing. The unique reconstruction conditions 
have been derived for normalized and unnormalized Laplacians 18,20 . As the held of 


graph signal processing is rapidly developing, we summarize some recent related works as 
follows. A least square approach has been proposed in jF] to reconstruct bandlimited graph 
signal from signal values observed on sampled vertices, using a centralized algorithm. An 
iterative method of bandlimited graph signal reconstruction has been proposed in |6], with 
the practical consideration of balancing a tradeoff between smoothness and data-htting. 
Two more efficient iterative reconstruction methods using the local set have been considered 
A necessary and sufficient condition for perfect reconstruction of bandlimited 


in 21,22 


graph signal has been derived in [T]. Readers may refer to Section 2.3 for more details of 
related works. 

Signal processing on graph is naturally related to distributed systems. For large-scale 
systems in lack of a central controller, e.g., sensor networks, distributed estimation and 
tracking 123 is an important topic. Algorithmic frameworks for distributed regression 124j 


and inference 25 have been studied to fit global functions based on local measurements in 


sensor networks. Consensus-based methods have been proposed in 26,27| to distributively 


compute the maximum likelihood estimate of unknown parameters. Distributed Kalman fil¬ 
tering has been introduced in |28| for target tracking of sensor networks. Diffusion RLS 129 


and LMS 30 algorithms have been proposed for distributed estimation over adaptive net- 
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works. To the best knowledge of the authors, there have been few works on the distributed 
reconstruction problem of bandlimited graph signal reconstruction. A related work [31 pro¬ 
poses an approximation method that calculates the graph Fourier multipliers distributively, 
which we will discuss in subsequent sections. 

In many practical distributed systems, a central processor is lacking and the majority of 
nodes have severely limited data processing power. Moreover, the node-to-node transmis¬ 
sion delay is non-negligible in large-scale networks, i.e., a given node cannot communicate 
globally with all other nodes to obtain instant fresh data. These difficulties with large 
distributed systems pose a new and practical challenge to graph-based signal processing: 
how to reconstruct graph signals distributively and efficiently? This is the motivation of 
the present paper, in which we propose a distributed algorithmic solution and answer the 
prior question to a reasonable extent. 

1.3 Contributions 

In this paper, we focus on the distributed recovery problem of graph-based time-varying 
signal. We propose a distributed algorithm, namely the distributed least square reconstruc¬ 
tion (DLSR), to adaptively reconstruct the missing values of a graph signal by allowing 
neighboring nodes to communicate with one another and make local updates. Due to the 
transmission delay caused by node-to-node communication, some out-of-band energy in¬ 
evitably occurs during the distributed signal reconstruction process. The DLSR algorithm 
uses a decay factor to dampen this out-of-band energy and achieves perfect reconstruction 
of the unknown graph signal. Theoretical convergence proof and error bounds of DLSR is 
given in our analysis. 

The rest of this paper is organized as follows. In Section II, some preliminaries are 
introduced, including the basics of graph signal processing and a review of existing related 
works. In Section III, the distributed reconstruction algorithm (DLSR) is proposed and 
described in detail. In Section IV, the proof of convergence and error bound analysis for 
the proposed algorithm is presented. In Section V, DLSR is evaluated using numerical 
experiments with synthetic as well as real-world data. 

2 Preliminaries 

2.1 Laplacian-Based Graph Signal Processing 

The concept of graph Laplacian is widely adopted in spectral graph theory [321 and signal 
processing on graphs [1|. For an undirected graph, the graph Laplacian (or unnormalized 
Laplacian, combinatorial Laplacian) is defined as 

L = D - A, 
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where A is the adjacency matrix and D is the diagonal degree matrix, whose elements are 
the degrees of the corresponding vertices. The normalized Laplacian is defined as 

C = D 


Both unnormalized Laplacian and normalized Laplacian are real symmetric positive-semidefinite 
matrices and all the eigenvalues are nonnegative 32 . In the rest of the paper, we mainly 


focus on the normalized Laplacian. However, we note that analogous results can be easily 
obtained for unnormalized Laplacian. 

In view of signal processing on graphs, the eigenvalues {Xk} of the Laplacian are regarded 
as frequencies and the corresponding eigenvectors {u/J are regarded as basis vectors. Con¬ 
sider an arbitrary graph signal f £ 9?^. Its frequency component corresponding to A k is the 
inner product between f and the eigenvector u&, denoted as 

N 


f {Xk) = (f,u fc ) = 


i— 1 

The eigenvectors associated with small eigenvalues have similar values on neighboring 
vertices, while the eigenvectors associated with large eigenvalues are the opposite. As a 
result, the frequency components associated with small and large eigenvalues correspond to 
the low-frequency and high-frequency parts of the signal, respectively 00 


Suppose f £ is a graph signal on an JV-vertex graph Q(V,£). We say f is u- 
bandlimited if its frequency components corresponding to eigenvalues larger than u are all 
zero. In other words, the spectral support of f is a subset of [0, w]. The subspace consisting 
of all cu-bandlimited signals on graph Q is called the Paley-Wiener space, which is a Hilbert 
space and denoted as PW U (Q). 

Suppose that for f G PW UJ {Q) only the entries on a selected set of nodes {/(it)}ues are 
known, where S C V is the sampled vertex set. The sampling and reconstruction problem 
is to recover the o;-bandlimited original signal f based on the sampled data {f{u)} u£ s- 


2.2 Frame and Signal Reconstruction 

Bandlimited signal reconstruction is closely related to the frame theory. We briefly introduce 
its basics in the following. 

Definition 1 A sequence of vectors {fjjiex is a frame in a Hilbert space ft, if there exist 
two constants 0 < A < B such that 

A\\t\\ 2 < ^|(f,f*)| 2 < H||f|| 2 , Vf eft. 

i£l 

Here the constants A and B are called frame bounds. 
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Definition 2 For a frame {f i}i£x, the frame operator S : % —> P is 

i£l 

where S is invertible and satisfies AI A S A BI. 

If the Euclidean matrix norm satisfies 111 — AS 11 < 1, then 

OO 

f = S _1 Sf = A ^(1 - AS) J 'Sf. 

j =o 

Consider the problem of reconstruction of an unknown signal f*. By defining 

k 

f( fc ) = A^(I-AS) J Sf*, 

3=0 

we can use the following iteration to iteratively reconstruct f* (see |34|), 

f (fc+i) = A Sf* + (j _ AS)f( fc ) = + AS(f* - (1) 

which achieves an exponentially shrinking error bound 

||f (A;) - f*|| < ||I - ASf ||f ( °) - f*||, Vfc > 0. 

2.3 Previous Works on Bandlimited Graph Signal Reconstruction 

We review some basic concepts and important results regarding band limited graph signal 
reconstruction, which have been established in existing works. 

Definition 3 f7gj / A set of vertices U C V(f?) is a uniqueness set for space PW i0 {Q) if 
Vf £ PW U (Q), f is uniquely determined by its values onU, i.e., Vf , g £ PW U (Q), f| u = g| u 
implies f = g, where f| u G ^ |W| is the restriction of f to the subset U. 

In order to perfectly reconstruct the bandlimited signals, we need the following relation 
between the sampling set and frame has been established. 

Theorem 1 fWj If the sampling set S is a uniqueness set for PW^Q), then {Vu,(S u )} u& s 
forms a frame in PW U1 (Q), and the upper bound B = 1, where Vuj{ ) is the orthogonal 
projection onto PW UJ {Q), and S u £ 3?^ is a graph signal on Q whose entries satisfy 

{ 1 , v = u: 

0 , v / u. 

A method called iterative least square reconstruction (ILSR) has been proposed to 
reconstruct the bandlimited signal iteratively. To be consistent with the results of our 
work, the method is rewritten as follows. 
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Theorem 2 |6j/ If the sampling set S is a uniqueness set for PW U {Q), then the original 

signal f* E PW U (Q) can be reconstructed using the sampled data {f*(u)} u£ s by the following 
ILSR method, 

f (fc+1) = +Vu (£(/*(“) - f {k \n))S u j , (2) 

Ws J 

where f ^ is a temporary result in the kth iteration. 


ILSR is derived from the method of projections onto convex sets (POCS) |35,36 , which 
is also known as the alternating projection method. The iteration ([2]) can be obtained by 
projecting onto the following two sets alternately, 


Ci ={f E lR N \f(u) = /*(u),Vu E S}, 
C 2 =PW U {Q). 


An equivalent derivation of ILSR can be obtained with the help of frame theory. Because 
of Theorem [I] the above method can also be obtained for A = 1 in the iteration 0 - 
Therefore, in 0 is the same as that in 0. 

In addition to ILSR, two more efficient graph signal reconstruction methods, namely 
the iterative propagating reconstruction (IPR) and the iterative weighting reconstruction 


(IWR), have been proposed and proved to be convergent 21,22 


Sampling and reconstruction of bandlimited graph signal is closely related to irregular 
sampling 37-40 or non-uniform sampling 1411 in the time domain, which sheds light on 
the analysis of graph signal. In fact, ILSR, IPR and IWR all have correspondence in time- 
domain irregular sampling, which are known as the Marvasti method 42 , the Voronoi 


method |43 and the adaptive weights method 37 


2.4 Notation 

For a graph Q and a cutoff frequency oj, PW U {Q) denotes the cj-bandlimited space of graph 
signal on Q . For any graph signal f E 9?^, V u ( f) denotes its orthogonal projection onto 
PW'jj(G), and TL,_|_(f) denotes the projection onto the orthogonal complement space of 
PW U (Q). The sampling vertex set is denoted as S, and the communication time delay 
between vertices u and v is denoted as t(u , v). It is assumed that one iteration is conducted 
at each time step. The true graph signal at the kth time step or iteration is denoted as 
fi fe \ whose entry associated with vertex u is fi k \u). ff^ denotes a biased estimate of fi^. 
In the reconstruction algorithm, f^ denotes the temporary result obtained after the /cth 
iteration, with f^ k \u) as its entry associated with vertex u. 
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Figure 1: An example of graph signal reconstruction over wireless sensor network. 


3 Distributed Reconstruction of Time-varying Bandlimited 
Graph Signal 

3.1 Motivation 


We consider the distributed reconstruction of a time-varying low-frequency signal defined 
over graph by sampling at a small portion of nodes. The problem could be described in the 
scenario of wireless sensor network (WSN). For a given WSN, there are unknown function 
f* \v) associating with node v at time t. Suppose that the function is slowly varying over 
both the time domain and the space domain. A snapshot of such function could be modeled 
as an unknown time-varying low-frequency signal fi^ £ PW U (Q ) located over a graphd] 
Suppose that the WSN is a hybrid network and only a small subset of nodes in S are 
equipped with sensors. As a result, one can measure the signal entries on support S as 


{/i^(^)}«e5 a t time t. Our purpose is to distributivedly estimate the function values at 
all nodes by using historical measurements at selected nodes {fi T \u)} u& s,r<t- 

See Fig. [I] for a demonstration of the raised problem. 

When the signal is time-invariant, the problem reduces to bandlimited graph signal re¬ 
construction, where the nodes with and without sensors correspond to sampled and missing 
data. In the distributed setting, the centralized iteration ([2]) no longer applies, as it is 
impossible for every node to obtain the instant estimation errors {f*(u) — /^( u )}«e.s- 

In what follows, we focus on the generalization of ILSR method to distributed systems 
and time-varying signals. We proposed an algorithm called distributed least square re¬ 
construction (DLSR). By letting each node conducting the iteration locally at each time 


1 One may argue how one could model the WSN as a graph, e.g., how the weighted edges are yielded 
among all nodes. However, this problem is beyond the topic of this paper. We always assume the graph is 
available or could be estimated by some available methods. 
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instant, DLSR can adaptively reconstruct the missing entries of a slowly time-varying graph 
signal. 

3.2 Algorithm Description 

The basic idea of DLSR is to spread the current estimation errors associated with the 
representative nodes (which are equipped with sensors) to all other nodes over the connected 
network. Every node iteratively updates its own estimation based on its received messages. 

The driver of the proposed algorithm is on those nodes with sensors, which calculates 
the error between the measurement f* k \u ) and the temporary estimation f( k \u) in the 
kth iteration by 

e (k \u) = fi k \u) - / (fc) (u), Vu € S. 

Then the estimation errors at node u (u E S) are transmitted to other nodes in the network. 
At the A;tli iteration, an arbitrary v collects a set of delayed but most recent estimation 
errors, 

{e R-rKD) (u)}ue5) WeV(S), 

where t(u, v) denotes the transmission delay from node u to node v |^} We denote the 
maximal transmission delay of the network by 

r= max t(u,v). 
ueS,veV{Q) 

Utilizing the most recent estimation errors, node v updates its local estimate by 

f {k+1) (v) = (1 - ^ k+l h +1 )f^ k \v) + Ii k+1 £ e ( fc -' r (“>U)( u ) (VojSu)( v) , W 6 V(g) (3) 

u£S 

where n k +1 and Pk+i denote the stepsize and decay factor, respectively. (VajS u )(v) denotes 
the entry at v of the lowpass component of S u , which could be calculated and stored before 
the system starts. 

Please refer to Table [l] and Table [2j which describe the detailed iterative process of 
signal reconstruction at the presentative nodes and the remaining nodes, respectively. 

3.3 An Example 

In order to provide some intuition for our distributed algorithm, let us refer to Fig. [l] and 
describe what happens on a typical node in the network. 

• As a representative node equipped with sensor, node 1 will get a measurement f* k \ 1) 
at the kth iteration and then calculate the estimation error e^Rl). The estimation 
error will be send to its neighbors of node 2, 4, and 5, and then forwarded to others. 
At the same slot, node 1 will receive the estimation errors of other nodes with sensors 


2 For simplicity, we may set t( - u,v ^(u) = 0 if k — t(u,v) < 0. 



Table 1: DLSR Algorithm at Representative Node u £ S. 


Parameter: S, {t{u' , u)} u ^ ueS , Hk, 

Initialization: f^°\u) = 0, calculate (V u d u /)(u),Mu' € S; 

For k = 0,1,2,--- 

1) Input: f* k \u)- 

2) Estimation: 

e ( fc )( w ) = fi k \u) - f (k \u); 

3) Communication: 

Send e^ k \u) and e ( fc - 1 -' r ( u '> u ))( lt / ) to neighbors, Mu' G S\u; 
Receive e^ k ~ T ^ u ' ,u ^(u') from neighbors, Mu' G S\u\ 

4) Update Storage: 

Save e (fc-r(«',«))( u /) ) v«' g 5\u; 

5) Update Estimation: 

/ (fe+1) («) = (1 ~ Mfe+i/5fe+i)/ (fc) (^) 

+mi £ («')(«)■ 

u'es 

End 


and use the most recent ones (e^ k 1 ^(2) from node 2 and e^ k 2 \3) from node 5). 
Consequently, it could update the estimation by 

/ (fc+ 1 )(l) = (l_ /Xjt+lA+ 1 )/W(l) 

+ Mfc+i ^e^(l) • (TU^i)(l) + e^- 1 ^(2)(P tJ ^ 2 )(l) + e^' _2 ^(3)(P w ^3)(l)^ • 

One may notice that the new estimation is not the output of the proposed algorithm, 
because our purpose is to estimate the strength of the signal associated with the node 
without sensor. However, the new estimation errors at representative nodes will be 
transmitted over the network to help all others to conduct their estimation. 

• As a regular node that is not equipped with sensor, at the kth. iteration, node 5 will 
receive e^ fc-1 ^(l), e^ k ~ 2 \2), and e( fe_1 )(3) from its neighbors. Then it will transmit the 
most recent estimation errors of nodes with sensors to its neighbors 1,3, and 6. The 
estimate of node 5 is updated by 

f^ +1 \5) = (l-fi k+1 p k+1 )f^(5) 

+ Mfe+1 ( e M(l)(P^ 1 )(5) + e( fc - 2 )(2)(P a; <5 2 )(5) + e( fe - 3 )(3)(P w 53)(5)). 
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Table 2: DLSR Algorithm at Non-representative Node v G V(£)\«S. 


Parameter: 5, {t(u, v)} ueS ,veV(G)\S, rtfc, Pk', 

Initialization: = 0, calculate {VuS u )(v),Vu G S ; 

For k = 0, 1,2, • • • 

1) Communication: 

Send e( fe - 1 - r ( u, ' y ))(u) to neighbors, Vu G 5; 

Receive e^ k ~ T ^ u,v ^ (u) from neighbors, Vu G 5; 

2) Update Storage: 

Save e^ k ~ T ^\u), Vu G 5; 

3) Update Estimation: 

/( fc + 1 )(u) = (l-/x fc+1 A +1 )/( fc )(u) 

+/Xfc +1 £ e^ k - T M\u)(VM(v); 

uGS 

4) Output: f( k+1 \v). 

End 


The new estimate will be sent out as a temporary result of the proposed algorithm. 

3.4 Discussions 

Asynchronization is one of the most common issues in distributed systems. In our dis¬ 
tributed reconstruction setting, asynchronization leads to communication delay between 
nodes that are connected through multiple links, which induces a deviation of the esti¬ 
mated signal from the bandlimited space. However, as long as the maximum delay in the 
network is bounded by a constant r, the proposed method can successfully annihilate the 
out-of-band estimation error and achieves perfect reconstruction. 

Node failure is also a common problem in WSN. The proposed DLSR is robust to both 
communication failure and sensor failure. In the former case, some links are broken and fail 
to work. The data packets have to be delivered through new route and the transmission 
delay may increase. Even in this case, the DLSR is still going to work, provided that the 
network remains connected and the maximum transmission delay is bounded. In the case 
of a sensor failure, some presentative nodes (i. e. u G 5) no longer obtain the sampled 
data, which means that they can act the same as the regular nodes. As long as the system 
is designed with some redundancy, the DLSR can still work, provided that there remain 
enough number of functional sensors. 

The proposed algorithm requires that the vectors {VojS u } u ^s be calculated in advance 
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and their entries be stored in the respective nodes. In some practical situation this pre¬ 
calculation could be unavailable. However, a distributed method proposed by [31 can be 
used to calculate {(7U<U)(u)} uei s approximately at node v. 

In fact, the operation TL(') for any given graph signal can be approximately calculated 
by a distributed method proposed by 31 . By this method, it will take some (depends on 
the network scale) rounds of data transmission and calculation to obtain the approximate 
projection. Therefore, another distributed method can be readily proposed by directly 
applying the above approximate projection to ILSR with a stepsize /r to track time-varying 
signal. Our method differs from this approach in the following aspects. 


• Supposing the period of a time step composed of one transmission and calculation 
is fixed in both methods, it will take some time steps to implement one iteration of 
ILSR by directly applying the method in [311 to ILSR. Therefore, the data used in the 
calculation are all sampled several time steps earlier. In our method, the iteration is 
conducted in one time step and use the current data at every node, which means the 
samples are as fresh as possible, the delay is only caused by transmission, and there 
is no waiting. Although nonuniform delays may violate the bandlimited property, 
we will prove in the following section that the out-of-band energy can be eliminated 
eventually. 


• The projection TL(') should be conducted for different graph signals as the itera¬ 
tion goes on by directly applying the method in [31 to ILSR, and each projection 
takes some time steps. In the proposed method, only pre-calculating (precisely or 
approximately using the method in |3l]) the frame elements, {Vu,S u } u& St are enough 
to reconstruct the graph signal, which is more economical. 


4 Convergence Analysis 

We will first study the convergence behavior of DLSR in a general situation that the ban¬ 
dlimited graph signal to be constructed varies slowly by time, and then specialize the result 
to a time-invariant case. In order to simplify the expression, we will fix stepsizes /i and (3 to 
be constants in studying the time-varying case. Finally, we let the stepsizes be diminishing 
and show that DLSR achieves a perfect reconstruction of time-invariant signal. 


4.1 Tracking Time-Varying Signal Using Constant Parameters 

The proposed DLSR algorithm is equivalent to the following iteration in the vector form 

f (fe+1) = (1 - M/3)f (fc) + M E ( F *« " F - fe) ) V ^ ( 4 ) 

u£S 


where 


4D 


- FW = diag - /( fc - T (“’*))(u)) 

l )i=l,-,N 
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is a diagonal matrix composed of the delayed estimation error at node u. 

Although V U S U is bandlimited for any u, by introducing , the delayed signal [Flu ~ 
Fu ' l )V IJ j8 u no longer belongs to the low-frequency subspace PW U] {Q). As a result, the 
sequence of estimated signals {f} are no longer cu-bandlimited. This existence of out-of- 
band energy makes DLSR critically different from its centralized version ([ 2 ]) and substan¬ 
tially complicates the convergence analysis. Since is not cu-bandlimited, we need to 
study its low-frequency and high-frequency components separately. 

For given node set S and cutoff frequency w, we may define an operator T on a graph 
signal f as 


Tf = v AY /w* 

\ues 

= Y fW'PuK- 


(5) 


u£S 


According to Theorem[lJ if S is the uniqueness set of graph Q with respect to w, {Vu8 u } u ^s 
is a frame in PW U (Q). For any f e PW UJ (Q), using the fact that 

/(«) = (V u f, 8 U ) = (f, V u 5 v ) f Vu e 5, 


Tf can be rewritten as 

Tf = E <f ,V u Su)V u S u , Vf ePW u (Q), 

u£S 

which is the frame operator of {V c j8 u } U £Si and the frame bounds are A and B. 
The definition of T implies 


I Tf 11 < 


and one has ||T|| < 1. By defining ff^ as 


Y 

ueS 


< llfll 


one further gets 


f* w = (/n + Tj^Tfr, 


Tfi fc) = /3fi fc) + Tfi fc) . 


( 6 ) 


(7) 


According to both Tfi^ and Tf**^ are within the low-frequency space PW U (Q). There¬ 
fore one can obtain from (7) that fi^ G PW w (Q). As a consequence, TE + fi^ = 0, where 
V UJ+ denotes the projection operator onto the high-frequency subspace which is the orthog¬ 
onal complement of PW U (Q). 

By defining the in-band error and out-of-band error as, respectively, 


e« = 


(k) 

e_i_ = 


7Ef (fe) - vJi k) 


V f ( fc ) _ V f ^ — V f 

/ oj+ 1 I o;+ A * — / uj+ 1 


( 8 ) 

(9) 
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the following proposition gives inequalities that and {e^} satisfy. Further, it will be 

shown that if the signal varies slowly enough, by properly selecting the stepsize p, DLSR 
can track time-varying signals. 


Proposition 1 Supposing the true signal satisfies 

fi k+1 \u) - f* k \u) < A, Vu € V(Q), k > 1, 
if A < min {A max ,/3B e+ /(\S\r)} and 


f^min ^ M ^ Hlin \ [1maxi 


Brn 


P + A 'C + \S |§r 2 A r 


the errors and {e^} satisfy 

e (k+i) < _ ^) e W + + M(/i)A, 

e( fe+1 ) < (1 - /i/3 - M)e (fc) + /i||T||ef } + ^ 2 C + M(p) A. 
/n f/ie above inequalities, 

M(p) = |5|^r 2 ^ 2 + \S\th + VN, 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 


A max is the positive root of 

|5|ir 2 (^4 m - |5|3^ A 2 + (2\S\rfiB e+ + A — (3 2 B 2 + = 0, 

Pmm and p max are the roots of 

(C + |«S|§t 2 a) p 2 + (|5|rA - fiB e+ ) p + VNA = 0, 

A and B are the frame bounds of T in PW U {Q), ||T|| is the norm of T, B v , B e and B e+ 
are constants satisfying 

(fi + A)B e = (f3+\\T\\)B e+ , (15) 

and C is a constant 


C = ((/? + ||T||)(R e + B e+ ) + B v ) . 

The proof of Proposition [I] is postponed to 7.1 


(16) 


Remark 1 According to Proposition [7J the out-of-band error (12) shows the necessity of 
the decay factor f3. If (3 = 0, the out-of-band error cannot be proved to converge, and the 
error may accumulate with the iterations. The decay factor (3 enhances the robustness of 
iteration. 

I (M 

Remark 2 The inequality (IS) shows that the out-of-band errore\ also affects the in-band 


error e^ k+1 \ which implies that the in-band error cannot be very small if the out-of-band 
error exists. In other words, it is important to eliminate the out-of-band error. 
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Taking the limit superior of (12) and (13), one obtains 


lim sup + — ) /z + 

k—yoo 


lim sup e^ < ( 1 + 

k—> oo 


E 
| T | 

T 


M(fi) 

Pp 

D/3 + E 


P + 


M(li) 


P + A iP + A)ix 


A 


(17) 

(18) 


where D and E are constants. The above results imply that for a constant stepsize /i, the 
(k) 

out-of-band error e\ will eventually get below a threshold that is determined by p and 
A. Similarly, will converge to the neighborhood of Vjt* k \ and the error is also 

controlled by p and A. 


Remark 3 Because bias is introduced by the multiple 1 — /z/3 in the iteration, the low- 
frequency and high-frequency components of the temporary estimate will get into the neigh- 

~(k) 

borhoods of ' and 0, respectively. Because of the influence of the decay factor (3, the 
reconstructed signal is biased. It will be proved in Corollary [I] that in the time-invariant 

~(k) 

case, these two components will exactly converge to V u f* ’ and 0. 


4.2 Recovering Time-invariant Signal Using Constant Parameters 

For the time-invariant case, fcan be written as f* and F*^ becomes f*(u) Ijv- Similar to 
(JgJ) , f* can also be defined as 

f* = ( / 9I + T)“ 1 Tf*. (19) 

Then Proposition [l] becomes the following corollary for A = 0. 

Corollary 1 For time-invariant true signal f*, supposing the in-band error and out-of-band 
error are defined as, respectively, 


e (k) = W'PJW 

ef = ||7U + f (fc) - V u+ l\\ = ||n,+f (fc) ||, 

the error {e^} and {e^} satisfy 

ef +1) <(l - iiP)e^ + p 2 C, 

e (k+i) < (! _ ^ - ^i) e ( fc ) + p\\T\\e {k) + i?C, 

„ . f 1 B r , PB e+ 1 

//C mm J, 

where the constants are the same as those in Proposition [7J 

Besides, © and (fl8|) become, respectively, 


lim sup < — p = [ D + — ) /z, 


C 


k—^oo 


P 


E 


(k) Dp + E f 

lim sup e (K) < — -— 1 

fc^oo - p +A \ 


p 


+ 


p 


/i. 


( 20 ) 

( 21 ) 


( 22 ) 

(23) 
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(k) 

Remark 4 For the time-invariant case, the out-of-band error e\ J will get below a threshold 
that is proportional to the stepsize p. It means that the out-of-band energy will be almost 
eliminated eventually along with the iteration if p is small. will converge to the neigh¬ 

borhood ofV^U, and its radius is also proportional to p. Therefore, for diminishing stepsize 
Pk approaching 0, V u f^ and V u+ i^ will strictly converge to and 0 , respectively, for 
a sequence of properly chosen diminishing stepsize. 


The bias will be estimated next. Because f*, f* E PW U] {Q), the iteration will converge 
to f*. Considering the definition of f* in (19), the bias satisfies 


f* - f* = (/3I + T) -1 Tf* - f* 

= (/3I + T) _1 Tf* - (/3I + T) _1 (/3I + T)f* 
= —/3(/3I + T) _1 f*. 


For f* E PW U (G), according to the frame bounds of operator T, we have 


I + T)- 1 fJ| < 


P + A 1 


and then 


If _ f || < 

I.* ^ 


p 


p + A" *"■ 

Thus, the bias is determined by p and decreases with the decrease of p. 


(24) 


Combining (22), (23), and (24), the following proposition gives the limit superior of the 
total error, as a function of P and fi. 

Proposition 2 The total error for the time-invariant case satisfies 

limsup ||f (fc ) — f* || < FP + G ^ + Hp, + Jpp, 

k —^oo P 


where F, G, H, and J are constants. 


Proposition [2] can be easily proved by summing up (22), (23), and (24). 


4.3 Recovering Time-Invariant Signal Using Variable Parameters 

Finally we will back to the general situation of variable stepsize and decay factor. According 
to Proposition [2} by discarding the higher order of a diminishing stepsize pk, the best Pk 
satisfies Pk ~ 0(y/~pk) to minimize the total error. Accordingly, the total error bound can 
be controlled by adjusting pk and Pk- In other words, the reconstruction error can be made 
arbitrarily small: DRSL achieves perfect reconstruction. 

The following proposition gives the convergence analysis for a special choice of dimin¬ 
ishing stepsize and decay factor. 
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Proposition 3 For diminishing stepsize fj,f. = \i\j\fk and decay factor (3^ = j3\/\fk, the 
total error satisfies the following inequality, 

||f (fc) -f*|| < K/<fk, 

where I\ is a constant. It means that the total estimation error decreases on the rate of 
l/\/k and converges to zero eventually. 

The proof of Proposition [3] is postponed to 

5 Experiments 

Experiments are designed to confirm the theoretical analysis and test the performance of the 
proposed distributed algorithm. The graph is generated by 100 randomly located nodes and 
the edges are generated by 4-nearest neighbors of the nodes, and the weights are inversely 
proportional to the square of geometric distance. Among the 100 nodes, 20 of them are 
randomly selected as the sampling set. The cutoff frequency is chosen to guarantee that the 
sampling node set is a uniqueness set, which can be determined by the method given in |[5]. 
[^] The bandlimited signal is generated by filtering the high-frequency components off. The 
transmission delay of each pair of nodes is simply regarded as the number of hops between 
them in the graph. The maximal transmission delay of this graph is 14. 

5.1 Tracking Time-Varying Signals 

5.1.1 Tracking Performance 

In this experiment, the tracking performance of DLSR is verified. The parameters are 
chosen as A = 0.005, /x = 0.1, and (3 = 10 -3 . The time-varying signal is generated by 
adding a random bandlimited increment whose largest absolute entry is A for each time 
step. The aiming signal and iterative results of four nodes are focused on, as illustrated in 
Fig. The nodes associated with the upper two subfigures are in the sampling set, and 
the nodes in the lower two subfigures are not in the sampling set. All the nodes can track 
the aiming signal for not dramatic changes. The proposed algorithm can track the slowly 
varying graph signal along with time. 

5.1.2 Parameters /3,/x, and A 

In this experiment, the regions for the parameters f3 and /x that guarantee the convergence 
of DLSR are plotted in Fig. [3] for different A, which describes the varying rate of time- 
varying signals. The experiment results show that for time-varying case the algorithm is 

2 Proposition 2 of [H]: The sampling set 5 is a unique set for PW^{Q) if the cutoff frequency satisfies 
uj < <7 m in, where cr 2 lin is the smallest singular value of (£ 2 ),sc, which is the submatrix of £ 2 containing only 
the rows and columns corresponding to the complementary set of S, and C is the normalized Laplacian of 

Q. 


7.2 
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Figure 2: Time-varying aiming signal and iterative results of four nodes. The proposed 
algorithm can track the aiming signal over time. 


not convergent if the stepsize is too large or too small. If the [i is too small, the estimation 
cannot track the varying signal. This will not happen for the time-invariant case (A = 0). 

5.2 Reconstruction of Time-Invariant Signal 
5.2.1 Convergence Performance 

In this experiment, DLSR is used to reconstruct time-invariant signals. The convergence 
curves of distributed and centralized algorithms with constant stepsizes /r = 0.01 and /i = 
0.02 are illustrated in Fig. [4| with the decay factor /3 = 0.01. The centralized algorithm 
uses fresh data from the sampled nodes, while the distributed algorithm uses data with 
transmission delay. It is easy to see that a larger stepsize results in a faster convergence. 
The bias caused by j3 can be seen in the convergence curves of distributed algorithms, while 
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Figure 3: The probability of convergence for different choices of /j and [5 in time-invariant 
and time-varying cases. 


the error of centralized algorithm shrinks exponentially to zero with no bias. 



Figure 4: The convergence curves of distributed and centralized algorithms with different 
constant stepsizes, where the decay factor is fixed /3 = 0.01. 
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Figure 5: The in-band and out-of-band errors for different (3 if the initial value has out-of- 
band energy. If [3 = 0 the out-of-band energy cannot be eliminated as the iteration goes, 
and it also leads to a relatively larger in-band error. For (3 > 0, a larger (3 will lead to a 
faster shrinkage of the out-of-band error, and a larger steady-state in-band error. 


5.2.2 In-Band and Out-of-Band Errors 

As proved in Corollary [I] both the in-band and out-of-band errors will decrease into a 
small bound as the iteration goes for f3 > 0. In this experiment, we set the initial value 
f (0 ) with about 10% of out-of-band energy and conduct the iteration with different decay 
factors (3 = 0,0.005,0.01,0.05, and 0.1. The stepsize is chosen as ji = 0.2. The in-band 
and out-of-band errors are illustrated in Fig. [5} The experiment result shows the necessity 
of the decay factor. Although there is no bias for (3 = 0, the out-of-band energy cannot 
be eliminated as the iteration goes. It should be noted that the curve for f3 = 0 is not 
comparable with the others because the out-of-band energy also affects the in-band error 
according to Remark [2| Since the out-of-band energy cannot be eliminated, it also leads to 
a relatively larger in-band error for j3 = 0. For (3 > 0, it can be seen from Fig. [5] that even 
though the initial value has out-of-band energy, it will shrink towards zero along with the 
iteration. A larger (3 will lead to a faster shrinkage of the out-of-band error, and a larger 
steady-state in-band error. 
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Figure 6: Convergence curves for different choices of (3 and /x . The decay factor (3 mainly 
determines the steady-state error and the stepsize [x mainly determines the convergence 
rate. 


5.2.3 Constant Parameters (3 and [x 

The convergence curves for different choices of constant (3 and /x are illustrated in Fig. 
[6j As analyzed above, the convergence rate is mainly determined by the stepsize /u. The 
steady-state error is composed of two parts, the bias and the steady-state error introduced 
by the constant stepsize. The latter is relatively small compared with the former, which 
is mainly determined by the decay factor (3. It is obvious that a smaller f3 will lead to a 
smaller steady-state error. 

The steady-state errors and convergence rates for different choices of (3 and /i are plotted 
in Fig. [7J For fixed (3, the steady-state error varies little with [x. It shows that the bias, which 
is determined by (3, is dominant in the total error, while /x influences the total error little. 
Since there is bias in the convergence, the convergence rate is approximately calculated as 
(||f( m ) _ f*||/1|f(°) — f*||) 1 / m , where m is the number of iterations when the error reaches 1.2 
times the steady-state error. The rate of convergence is smaller for larger /x, which means 
it converges faster. 

5.2.4 Diminishing Parameters /x^ and (3^ 

An experiment for diminishing stepsizes and decay factors is conducted and the convergence 
curves are shown in Fig. |ij The stepsizes are chosen as [Xk = [X\/\fk with /xi = 0.05 or 
0.02. The decay factor are f3k = j3i/y/k with (3\ = 0.1 or 0.01. All the curves decline along 
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Figure 7: The steady-state errors and convergence rates for different choices of /3 and ^. 


the iteration. Among the four curves, it can be seen that a larger m and a smaller /3± may 
lead to a faster convergence. 



Figure 8: The convergence curves for diminishing stepsizes and decay factors in Proposition 

El 


5.3 Experiments with Real Data 


The sensor network data of Intel Berkeley Research Lab 44 is used in this experiment. The 
data is collected from 54 sensors in the lab and sampled every 30 seconds from February 
28th, 2004, including temperature, humidity, light, and voltage. In our experiment, the 
graph signal is composed of the temperature of the sensors. We extract the data from 
01:06:15 to 17:56:15 on February 28th, 2004, during which time there is less missing data. 
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Figure 9: The temperature of each node in the sensor network data and the relative error 
of DLSR. 


Taking time and space smoothness into consideration, the missing data is interpolated by 
conducting the MATLAB function scatteredlnterpolant with all the existing data. Then 
the completed data is regarded as the original time-varying graph signal. The graph is 
established by the 4-nearest neighbors of the positions of the sensors, and the weights are 
inversely proportional to the square of geometric distance. We randomly choose 20 sensors 
and reconstruct the temperature of the other sensors. By selecting /r = 0.1 and f3 = 10 -3 , 
the time-varying graph signal is reconstructed by DLSR. The temperature of each node and 
the relative error are illustrated in Fig. [9j The steady-state relative error is around 3%, 
which verifies the effectiveness of DLSR. 

6 Conclusion 

In this paper, the distributed least square reconstruction algorithm (DLSR) is proposed to 
estimate and track the unobserved data of a time-varying graph signal adaptively. The low- 
frequency and high-frequency components of the recovered signals are theoretically proved 
to converge, respectively, to their true values. The out-of-band energy caused by node-to- 
node transmission delay can be eliminated by using the decay factor, which introduces a 
controllable bias. The expression of the overall error bound is given as a function of the 
stepsize and decay factor, and can be made arbitrarily small. Numerical experiments on 
both synthetic and real world data verify the performance of the proposed algorithm and 
show that DLSR is able to track slowly varying graph signals adaptively. 
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7 Appendix 


7.1 The Proof of Proposition [l] 

First, besides the in-band error and out-of-band error defined in ([ 8 ]) and Q, two sequences 
of quantities are introduced as 


6 {k) 


f(fc) _ f(fc-i) 


(25) 


r,W 


u£S 


(26) 


Then we will prove the following inequalities, where the proofs are postponed to the end of 
this subsection. 


e (fc+i) <(i -up- fiA)eW + v\\T\\e {k) + W {k) + (VN + jti|«S|r) A, (27) 

e (k+i) _ ^y{k) + ^(fe) + ^ + /i |«S| T ^ A, (28) 

T—1 

7 (fe) <>/i5iE tf(fc_i) . (29) 

4 = 0 

<5 (fc) </r ((/3 + 11T11) + ^|5|rA, (30) 

Plugging (30) into (29), we have 

</i(c + |«S|§T 2 A) , (31) 


where C is dehned as (16). 

By plugging (31) into (28) and (27), respectively, we have © and ( pT3] ) readily. As a 
consequent, if we could demonstrate that {eW}, {e^}, and { 77 W} are bounded by some 
constants, respectively, the proof of Proposition [l] will be closed. In what follows, we will 
prove this by mathematical induction. 

For i = 1, these quantities are obviously bounded. We will then prove that if the 
preceding k — 1 items of {eW}, {e^}, and { 77 W} have respective bounds of B e , B e+ , and 
B v , their feth items are also bounded by the same limits. 


1. According to (31), rf k ^ is bounded by B v 


2. Plugging ^ < B e+ in (12), we have 


e {k) < (1 - pP)B e+ + \x 2 C + M(n) A, 


where M(/r) is dehned in 


We then have < B e+ if 


(c + |cS| 5 t 2 a) f! 2 + (|<S|rA - pB e+ ) /r + VNA < 0 . 


(32) 


23 



















One may notice that the inequality (32) has solutions for /x, if and only if the following 
inequality holds for A, 

(|«S|rA - (3 B e+ f - 4\/iVA (c + |S|ir 2 A) > 0, 
which is equivalent to 


|5|ir 2 (an* - |S|3) A 2 + (2\S\T/3B e+ + 41\^<?) A < (3 


2n2 

B e+ - 


(33) 


Because its left hand side is an increasing function of A, the inequality (33) is satisfied 


if A < A max , where A max can be solved from (33). 


Then from (32) the range of /r can be determined as /x m j n < /i < /U max , where both 


^ min and /i max are related to A. 

3. Plugging e < - fc ~ 1 ' ) < B e and e+ _1 ' ) < B e+ into (27), we have 


e (fc) < (1 -nfi- nA)B e + fi\\T\\B e+ + n 2 C + M(/i) A. 


Using (15), we can also obtain that e^ k > < B e if (32) is satisfied. 


Consequently, {e^}, {e^}, and {rj^} are bounded, and then Proposition |lj is proved. 


To end this subsection, we will prove (27), (28), (29), and (30). To simplify the ex¬ 


pression, we introduce two vectors to denote the misalignment of estimated signal and the 
increment of true signal by 


d (fc) = f(*0 _ fW ; 

n (k) _ p(k) _ p(k— 1 ) 

u - isk 1 


(34) 

(35) 


and two diagonal matrices to denote the errors at node u caused by delayed true signal and 
estimated signal by 


E& } = I N ~ F$, 
E ( fc) = fW{u)I N - F^, 


(36) 

(37) 


7.1.1 The Proof of 


and (28) 


According to @, ([5j), and Q , we have 

d (fc+1) =(1 - tf) d( fc ) - c( fc+1 ) - iipH k) + n ^(Fi fe J - FW)7>A 

u£S 

=(1 - /r/3)d( fc ) - c( fc+1 ) - /iTd^ + ^T(f< fc ) - fi fc) ) + /x ^(F^ - fW)PJ u 

u£S 

=Q ( V ,d« + P w+ d (fc) ) - c ( fc+1 ) + M Y, ( 38 ) 


u£S 


u£S 
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where 


Q = (1 — hP) 1 n — mT. 

Considering the definition of T in ([5]), for any f, we have 


VMVJ = QVJ, 

(39) 

VMVco+f = -liTP u+ f, 

(40) 

V u+ QV u f = 0, 

(41) 

Pto+QPu+f = (l-nP)V u+ f. 

(42) 


Therefore, according to (39) and (40), the low-frequency part of (38) is 
iP w d( fc+1 ) = QV^dW - pTP u+ dW - iP w c( fc+1 ) + ptV w Y ~ vPu Y E *«P<A- 

(43) 


itSiS 


u&S 


Since the frame bound of {VuPu}ueS satisfies A < B < 1, if we choose a stepsize [i 
satisfying (j, < 1/(P + A), the assumption AIn A T A BI n implies, note that V cu d < - k ' ) £ 

PWUG), 

IIQH < 1 — /r/3 — )jlA. 


According to (10), 


c( fe+1) 


< 


,(k+l) 


< Vna. 


The definition of Fimplies 


—rAIjv A A tAIat, 


and then the last term of (43) is bounded by 


V I E^V 8 

• lj I / J J - J *n / u u ' 


KuES 


< |<S|rA. 


Taking the norm of (43) and combining the above inequalities, the inequality (27) is 
obtained. 


According to (41) and (42), the high-frequency part of (38) is 

7V d(fe+1) =(1 - ^P)Vu + d^ - + nV u+ Y - fiV u+ Y E&V.A 


u£S 


u£S 


Following the samilar approach of proving (27), the inequality (28) is proved. 
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7.1.2 The Proof of (29) 


By the definition of F„ , 

K k) r.s u f =J_] - fV*- T M\u )| 2 \(v u s u ){v)\ 

V ^ 

< E iwa)«i 2 (e |/ ( ‘ -i) (“> - / , ‘ _i_1) (*‘)|') 

v V i=0 / 


where 


\i =0 


X ] l (^)(^)| 2 = ll ^ l | 2 < 1 - 


Therefore, 


v (k) < 


E || E i fe) ^ 


u£S 


T— 1 


uGtS 7=0 

Ev4slfE|/ ( ‘-”w-/ ( ‘- i - 1, ( 


< 

7=0 \u£S 

i =0 


which is (29). 


7.1.3 The Proof of 

According to Q, ([5]), and Q, 


f(fc) _ f(fe-i) = _ M (Fit _1) - F^- 1 )) V U 6 U 

u£S 

= - /r/Jd^" 1 ) - /xTd^ 1 ) + fi E Fit _1 VuA - H E 

u£cS m£cS 

= - M (/3Ijv + T)d( fc “ 1 ) + ,i E - !i E 


U£cS 


u£S 


Taking the norm of (44), the inequality (30) is obtained 


(44) 
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7.2 The Proof of Proposition [3] 

According to the proof of Proposition [lj setting A = 0 and using ||T|| < 1 and (3 k < /3i, the 
inequalities (27)-(30) for variant {fik } and {/3k} becomes 

e (fc) <(1 - HkPk - PaP 4)e (fc_1) + n k e+^ 1) + Pfc?7 (fc-1) 

ef <(1 - Vk/3k)e+~ 1) + /W fc_1) 
v {k) 


i =0 


A <fi k ((/3i + l)(e (fc - 1} + 4 fc - 1} ) + r^ 1 ) 


(45) 

(46) 

(47) 

(48) 


Similar to the proof of Proposition |TJ it is easy to see that {e^}, {e^} and {?/*)} are 
bounded by constants B e , B e+ and B v for = f-ii/Vk. Plugging (47) and (48) into (46), 
we have 

e+ } < (1 - 7tfc/3fc) e +~ 1) + C"pfcPfc-r, (49) 

where C' = Ty/\S\(2{B e +B e+ ) + B^}. For n k = fii/y/k and (3 k = /3\/\fk, if ^ < 
L+/y/k — 1 is satisfied for a constant L +1 according to (49), < L + /\fk as long as the 

following inequality is satisfied, 


( _ Pill \ L + 

V Vk</kJVk^T 


Pi 


The inequality above is equivalent to 
fijC' yfkyfk^l 


+ 


y/ky/k — T 

Vk 


< 


L± 

yfk 


L+ y/k — t /y/k + y/k — \){y/k + y/k — 1) 


< Pi/3i. 


(50) 


Because the first term of the left side approaches n\C'/L+ and the second term approaches 
0 when k is large enough, by selecting a constant L + appropriately, (50) is established and 
then we have < L + / yfk. 

Plugging (47) and ( |48| ) into (45), we have 

e {k) < (1 - HkPk ~ Pfc-4)e (/c_1) + Hk e +-^ + CWPfc-T- 


If e^ k 1 - 1 < L/y/k — 1 is satisfied for a constant L, e^ < L/yfk as long as the following 
inequality is satisfied, 

^ i a\\ L i ^ L+ i r> d < A 
V VkKVk )) y/k^l Vk y/k^l yfkyfk^r ~ yfk' 

which is equivalent to 

mL+ . fia #k =t _ n _ ( (h\ 

L + L ({/k+^k^T)(Vk + Vk^T)~ IH \ Vk)' 
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Both the second and third terms of (51) approach 0 when k is large enough. By selecting 
a constant L appropriately, (51) is established and then we have e^ < L/\fk. 


Based on the analysis above, we have 

||f (fc) -f*|| < 


Pk 


< 


Pk + A 

fh 

A 


f*|| + e, 

f* II + L+ + L 


( fc ) i Jk) 


1 

w 


(52) 


when k is large, and Proposition [3] is proved. 
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